# -*- coding: utf-8 -*-

import pandas as pd
import numpy as np
import seaborn as sn
import matplotlib.pyplot as plt
import datetime

# params
params = {
    'legend.fontsize': 'x-large',
    'figure.figsize':(30,10),
    'axes.labelsize': 'x-large',
    'axes.titlesize': 'x-large',
    'xtick.labelsize': 'x-large',
    'ytick.labelsize': 'x-large'
}
sn.set_style('whitegrid')
sn.set_context('talk')

plt.rcParams.update(params)
pd.options.display.max_colwidth = 600
from IPython.display import display, HTML

'''
# read data
train = pd.read_csv('bike/day.csv')

print(train.info())
print(train.describe())
# 类别型特征
categorical_features = ['season', 'yr', 'mnth', 'holiday', 'weathersit',
                        'weekday', 'workingday']
# categorical_features = ['season', 'mnth', 'weathersit', 'weekday']
for col in categorical_features:
    print ('\n%s属性的不同取值和出现次数'% col)
    print(train[col].value_counts())
    train[col] = train[col].astype('object')
# 数值型特征
numerical_features =['temp', 'atemp', 'hum', 'windspeed']
# a= train[numerical_features].hist()
# print (a)
# plt.show()
# 特征与目标之间的关系
# 骑行年
# sn.violinplot(data=train[['yr','cnt']],x='yr', y='cnt')
# plt.show()
# 一年中每天

train['date'] = pd.to_datetime(train['dteday'])
train['dayofyear']= train['date'].dt.dayofyear
# print(train.head)
# fig,ax = plt.subplots()
to_data = train[["dayofyear", "yr", "cnt"]]
print(to_data.head())
print(to_data.describe())
# plt.figure(figsize=(10,3))
sn.lmplot(x="dayofyear",
             y="cnt",
             data=to_data,
             hue ="yr"
          ,fit_reg=False
             # hue_order =['0','1'],
#              ax=ax
            )
# hue='yr',
plt.title("daily distribution of counts")
# 季节
sn.violinplot(data=train[['season','cnt']], x='season', y='cnt')
# plt.show()

fig, ax = plt.subplots()
sn.barplot(data=train[['season', 'cnt']], x='season', y='cnt')
ax.set(title="Seasonly distribution of counts")
# plt.show()
# 月份
sn.violinplot(data=train[['season','cnt']], x='season', y='cnt')
# plt.show()
plt.clf()
fig, ax = plt.subplots()
sn.barplot(data=train[['mnth', 'cnt']], x='mnth', y='cnt')
ax.set(title="Monthly distribution of counts")
# plt.show()
# 天气
plt.clf()
fig, ax = plt.subplots()
sn.barplot(data=train[['weathersit', 'cnt']], x='weathersit', y='cnt')
ax.set(title="weathersit distribution of counts")
# plt.show()
# 工作日节假日
plt.clf()
fig, (ax1, ax2) = plt.subplots(ncols=2)
sn.barplot(data=train, x='holiday', y='cnt', ax=ax1)
sn.barplot(data=train, x='workingday', y='cnt', ax=ax2)
# plt.show()
plt.close()
# 数值型特征与y的相关性
corrMat = train[['temp', 'atemp', 'hum'
                 ,'windspeed', 'casual','registered','cnt']].corr()
mask = np.array(corrMat)
mask[np.tril_indices_from(mask)] =False
sn.heatmap(corrMat, mask=mask, vmax=.8, square=True, annot=True)
plt.show()
'''

'''
# 特征工程
import pandas as pd
import numpy as np

train = pd.read_csv("bike/day.csv")
# print(train.head())
# for col in list(train):
#     print(train[col].value_counts())
# 对类别型特征，观察其取值范围及直方图
categorical_features = ['season', 'mnth', 'weathersit', 'weekday']

# 数据类型变为object，才能被get_dummies处理
for col in categorical_features:
    train[col] = train[col].astype('object')

X_train_cat = train[categorical_features]
X_train_cat = pd.get_dummies(X_train_cat)
# print(X_train_cat.head())

# 最大值最小值标准化
from sklearn.preprocessing import MinMaxScaler
mn_X = MinMaxScaler()
numerical_features =['temp', 'atemp', 'hum', 'windspeed']
temp = mn_X.fit_transform(train[numerical_features])
# print(temp)
X_train_num = pd.DataFrame(data=temp, columns=numerical_features, index=train.index)
# print(X_train_num.head())
# 合在一起
X_train = pd.concat([X_train_cat, X_train_num, train['holiday'], train['workingday']]
                    ,axis=1, ignore_index=False)
# print(X_train.head())
FE_train = pd.concat([train['instant'], X_train, train['yr'], train['cnt']]
                     ,axis=1)
print(FE_train.head())
print(FE_train.info())
FE_train.to_csv('bike_day_FE.csv',index=False)
'''

from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error  #评价回归预测模型的性能
from sklearn.linear_model import LinearRegression
from math import sqrt

# 读取去量纲后的数据
df = pd.read_csv("bike_day_FE.csv")
y = df["cnt"]
X = df.drop(["cnt", "instant"], axis=1)
feat_names = X.columns

X_train, X_test, y_train, y_test = train_test_split(X,y,random_state=2,test_size=0.2)
# print(X_train.shape)
# print(X_train.head)


#最小二乘线性回归
# 1. 线性回归
lr = LinearRegression()
lr.fit(X_train, y_train)
y_test_pred_lr = lr.predict(X_test)
y_train_pred_lr = lr.predict(X_train)
# 权重系数
# fs = pd.DataFrame({"columns":list(feat_names), "coef": list((lr.coef_.T))})
# print(fs.sort_values(by='coef', ascending=False))

def RMSE(targets, preds):
    return sqrt(mean_squared_error( targets , preds))
# RMSE 评估模型性能
print('The RMSE of LinearRegression on test is %.6f'  % RMSE(y_test, y_test_pred_lr))
print('The RMSE of LinearRegression on train is %.6f'  % RMSE(y_train, y_train_pred_lr))
# 观察残差分布，看是否符合：噪声为0均值的高斯分布
f,ax = plt.subplots(figsize=(7,5))
f.tight_layout()
ax.hist(y_train-y_train_pred_lr, bins=40, label='Residuals Linear', color='b', alpha=0.5)
ax.set_title("Histogram of Resicuals")
ax.legend(loc='best')
# plt.show()
# 观察预测值与真值的散点图
# plt.clf
plt.figure(figsize=(4,3))
plt.scatter(y_train, y_train_pred_lr)
plt.plot([0,1],[0,1], '--k')
plt.axis('tight')
plt.xlabel('true price')
plt.ylabel('predicted price')
plt.tight_layout()
# plt.show()


# 2.正则化的线性回归
from sklearn.linear_model import RidgeCV
# 超参数范围
alphas = [0.01, 0.1, 1,10,100]
# 生成ridge实例
ridge = RidgeCV(alphas=alphas, store_cv_values=True)
# 模型训练
ridge.fit(X_train, y_train)
# 预测
y_test_pred_ridge = ridge.predict(X_test)
y_train_pred_ridge = ridge.predict(X_train)
# r2_score 评估模型性能
print('The RMSE of RidgeCV on test is %.6f'  % RMSE(y_test, y_test_pred_ridge))
print('The RMSE of RidgeCV on train is %.6f'  % RMSE(y_train, y_train_pred_ridge))

# 可视化
# print(ridge.cv_values_)
mse_mean = np.mean(ridge.cv_values_, axis=0)
print(np.log10(ridge.alpha_)*np.ones(3))
plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas), 1))
plt.xlabel('log(alpha)')
plt.ylabel('mse')
# plt.show()
print('alpha is ', ridge.alpha_)

from sklearn.linear_model import LassoCV
alphas = [ 0.01, 0.1, 1, 10,100]

lasso = LassoCV(alphas=alphas)
lasso.fit( X_train, y_train)

lasso.fit(X_train, y_train)
y_test_pred_lasso = lasso.predict(X_test)
y_train_pred_lasso = lasso.predict(X_train)
# RMSE 评估模型性能
print('The RMSE of Lasso on test is %.6f'  % RMSE(y_test, y_test_pred_lasso))
print('The RMSE of Lasso on train is %.6f'  % RMSE(y_train, y_train_pred_lasso))
# 观察预测值与真值的散点图
plt.close()
plt.figure(figsize=(4,3))
plt.scatter(y_train, y_train_pred_lasso)
plt.plot([0,1],[0,1], '--k')
plt.axis('tight')
plt.xlabel('true price')
plt.ylabel('predicted price')
plt.tight_layout()
# plt.show()

mses = np.mean(lasso.mse_path_, axis=1)
plt.plot(np.log10(lasso.alphas_), mses)
plt.xlabel('log(alpha)')
plt.ylabel('mse')
# plt.show()







